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1 Introduction 

Taylor’s formula shows how to approximate a certain class of functions by 
polynomials. The approximations have two nice properties. They are ar- 
bitrarily good in some neighborhood whenever the function is analytic and 
they are easy to compute. Our goal is to give an efficient algorithm to ap- 
proximate a neighborhood of the configuration space of a dynamical system 
by a nilpotent, explicitly integrable dynamical system. For a class of dynam- 
ical systems analogous to the analytic functions, this approximation will be 
arbitrarily good in some fixed neighborhood and easy to compute. 

In [2], we give an algorithm which given a rank r yields two vector fields 
E\ and E2 on with the properties 

1. The vector fields E\ and E 2 generate a Lie algebra isomorphic to the 
free, nilpotent Lie algebra on 2 generators of rank r. Let n denote the 
dimension of this Lie algebra. 

2. If Ei, for i = denotes the vector fields corresponding to the 

Hall basis of a free nilpotent Lie algebra, i.e. £3 = £4 = 

[E 3) Ei], E s = [£3,^1, efc., then 

s(o)= ^? ,=i "■ 
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3. The trajectory t — ► y(t) satisfying 

y(t) = «i(f )Ei(y(t)) + u 2 (t)Eh (y(0), y(0) = y° € R n 

can be written explicitly in terms of quadratures involving the func- 
tions t — ► ti(t). 

Let £ denote the configuration space for this system. We sometimes refer 
to this system as the model system. 

Suppose we are given an arbitrary system of the form 

i(t) = m^F^xit)) + u 2 {t)F 2 (x{t)) i 

where x(t) € R fc , F\ and F 2 are vector fields on R fc , and t -*• u;(i) are given 
measurable controls. Let T denote the configuration space for this system. 
In order to approximate this control system to m t/l -order, we compute the 
iterated Lie brackets Fj corresponding to the Hall basis F,-, obtained by 
substituting F\ and F 2 for E\ and £* 2 - 

Our main construction is an approximating map A which maps the con- 
figuration space £ of the nilpotent model to the space T with the property 
that the images under A of trajectories in £ with measurable, bounded con- 
trols stay close to their counterparts in F, provided that the Fi satisfy a kind 
of analyticity. The exact analyticity we require is described at the beginning 
of §3. 

The map A turns out to be a polynomial map from R n to R*. It depends 
only upon the vector fields E{ and F* and not upon the particular controls 
ifj(t). For this reason, it can be precomputed and used to compute efficiently 
a tubular neighborhood of trajectories around a given reference trajectory 
x(t). 

Nilpotent Lie algebras have been an important tool in control systems 
beginning with the work of Krener and Hermes, see [7], [4], [5], and [1]. The 
point of view of these papers was to use nilpotent Lie algebras in order to 
obtain theoretical results about properties of control systems. The point 
of view here is to focus on the some of the computational aspects of using 
nilpotent normal forms. In particular, we give an efficient algorithm to 
compute the map A mapping trajectories of a model nilpotent system to 
a given system and an algorithm to integrate a tubular neighborhood of 
trajectories around a given fixed reference trajectory. 

Section 2 defines the approximating map A. Section 3 defines the an- 
alyticity required for our algorithm. Sections 4 and 5 state and prove the 
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main theorem. Section 6 shows how to apply the main theorem to integrate 
simultaneously many trajectories around a given fixed reference trajectory. 
The final section contains some examples. The appendix contains the Math- 
ematica code we used to compute the examples. 

For unexplained terminology involving Lie algebras and Hall bases, see 
[3] and [6]. 

2 An approximating map 

We begin with an informal description of the main idea. The vector fields Ei 
define a map <$£;:R n — ► R n by (si,...,s n ) — ► exp(£r=i «*£*)• This takes 
the tangent space T(R n ) and flows it out into R". — * R* is defined 

similarly, by exponentiating £1 s t *F^ The lambda map is then ^f 0 ^ 1 The 
map is always invertible, so A is always defined. But since $f may be 
non- invertible, for instance, if k < n, then A may also be non-invertible. 

Computing <Pp is usually just as hard as solving for an arbitrary trajec- 
tory, so we take approximations of $f instead. Since our approximation is 
supposed to be good to order m, we use m applications of a Picard iteration 
scheme. To 0 th order, x° is 0. Substituting this into i 1 (r) = £$iFi(0) we 
get the first order approximation x*(r) = Next, we would like to 

solve x 2 (r) = £ SiFi(x l (r)). This is not usually solvable explicitly, but as 
we’re only interested in a second order approximation to the solution, we 
can take a first order power series expansion in the flow parameter r. Thus, 
we solve x 2 (r) = a 0 (t) + ai(f)r. We repeat this process and, eventually, get 
x m (l) to agree with $»f to order m. 

More precisely, to state and prove our main theorem requires the follow- 
ing four definitions. 

Definition 1 For fixed S\ , . . . , $ n , write x(r; , . . . , s n ) for the solution x(r) 
of x(r) = 53«tFf(«( r ))» *(0) = 0, and define operators by 

¥(z(r;si,...,s n )) = f ^2siFi(z(r] Sl ,...,s n ))dr 

Jo • 

W'(z(r; «!,...,<»)) = / T j (52 s i F i(x(r; *i, •••,*„))) dr, 

Jo i 

where T J represents the j th -order Taylor approximation with respect to the 
variable r . 


3 


Definition 2 The trajectories x j (r; s x , . . . , s n ) are defined inductively: 

a;0 ( r i S 1 1 »*n) = o, S 1, . . • , s n ) = *1 1 I s n))- 

Definition 3 The map $£(si, . . . ,s„): R n -» R "»'« defined as the time one 
flow of the trajectory x m , namely x m (l; «i , • • • , s n ). 

Definition 4 The m** 1 -order approximation to the A map, A m , is defined by 

A m = <f>J? o . 

3 The generalized Baker-Campbell-Hausdorff formula 

Proving convergence of the algorithm, requires generalizing the Baker-Camp- 
bell-Hausdorff formula (BCH). The BCH formula writes the product of two 
exponentials (the composition of two flows) as the exponential of a series 
in the brackets (a constant flow involving the higher order brackets of the 
vector fields). A trajectory in a control system is a limit of compositions of 
piecewise constant flows, and we can use the BCH to derive a constant flow 
involving a series in the higher brackets which arrives at the same point. At 
the formal level and for systems whose vector fields generate nilpotent Lie 
algebras, BCH holds exactly, providing computable “geodesic normal coor- 
dinates”. To approximate trajectories in other control systems, we must 
assume that BCH converges for them as well. We give an example later 
where BCH does not converge, but the vector fields are not analytic. To our 
knowledge, no general criteria are known which imply convergence of BCH. 

We now introduce the analyticity requirement we need in order to prove 
convergence of our algorithm. We say that two vector fields are BCH analytic 
in case there is a 6, with 0 < 6 < 1, such that any Lie algebra elements X\ 
and Xi in the Lie algebra generated by the F; of the form |Ai| < |Aj| < 6, 
satisfy the following estimate 

|e**e*» - e c ‘<*>.**)+ ..•+*»(*,, X,)| < tt ”*|A 1 ||X 2 r, (1) 

where c,(X i,X 2 ) are all the terms of weight i in the BCH formal expan- 
sion and a is a positive constant depending on the Lie algebra. The term 
a m |Xi ||X 2 | m comes from combining the fact that a converging series is geo- 
metric in a smaller ball, and hence a bound on the next term is comparable 
to the error, and the fact that every non-vanishing Lie element of weight 
> m has at least one Aj (the smaller) factor in every term. 
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Varadarajan [8] shows that the estimate |[X, Y]| < M |X||Y[, for all X, Y 
in the Lie algebra generated by F\ and F2, implies Condition (1). Note that 
this estimate holds if F\ and F 2 generate a finite dimensional Lie algebra. 

Lemma 5 Given a positive integer m f let n be the dimension of the free 
nilpotent Lie algebra of rank m on two generators. Given u,-(t) and [0 , T] 
as abo ve f then there exist constants si,...,s n sn ch given vector fields 
F\ and F 2 as above , and the higher brackets F3, ... } F n corresponding to the 
Hall basis elements , then the trajectories 

x(t) = «i(t)Fi(x) + u 2 (t)F 2 (x) t x(0) = 0 


and 

z(t) — s\Fi(z) + . . . + s n F n (z), y( 0) = 0 

satisfy \x{T)-y(T)\ < (aT) m+I . 

Note that the s* depend only on the ti» and not on the F{. The case tii(t) = 
u 2 (^) = T — 1 is the standard BCH formula. 

Proof. Let v\(t) and ^(f) be step-function approximations to u\(t) and 
ti2(0> respectively, chosen well enough so that the trajectory with the t>j(t) 
as controls stays t m+1 close to x(t) over the interval [0,T]. Suppose that the 
Vi(t) are constant except at the times (t\ , . . . , t/g). To start, consider the flow 
z l (t) with controls v t -(t) on the time interval [0, til- Trivially, there is a flow 
w l (t) with constant coefficients ]jT£= x sjF, such that w l (ti) = namely 

s] — Vi . The flow z 2 (t) is the flow u? 1 ^) for t £ [0,ti] followed by the flow 
with controls Vi(t) for t £ (ti,ts]. By the BCH convergence condition (1), 
there is a flow w 2 (t) with different constant coefficients s 2 Fi such that 
\w 2 (ti)-z 2 (ti)\ < (t 2 — ti)(aT) m . At stage j, the error is < (tj+i — t ; j(aT) m , 
and so the total error is < (aT) m — tj) = (aX , ) m+1 as desired. ■ 

Since the BCH formula holds exactly in the configuration space E, the 
trajectory z(t) is exactly A(y(f)) in this case. 

4 The Picard-Taylor Method 

Recall that the map $jp is the easily computed Picard-Taylor approximation 
of the exponential map 4>/r(s) = 5Z s%Fi. In this section, we prove that the 
error introduced by using is 0(s m+1 ), just as in Taylor’s theorem. 

The convergence of a Picard iteration scheme depends on a bound B and 
a Lipschitz constant L for the vector fields: |F,(p)| < B and |F,(p) — F,(g)| < 
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L\p — g| for p and q in some ball about the origin. The choice of ball 
determines B , L and a fixed time T' , such that the approximate trajectories 
do not leave the ball where these estimates hold. 

Lemma 6 Given Fi as above, there is a constant 6? > 0 and a positive 
constant C such that |(SJ? - $> F )(s)| < Cjs | m+1 «>Aenever |s| < 62 - 

Proof. Let z°(r ,s„)) = 0 and z>+ l (r,s u . . . ,s n )) = Exam- 
ine the convergence of z J (l;«i, . . . , s n ) to $f(*) = s *^») as j -+ 00 . 

Assume, as an induction hypothesis, that z,(r;s) and Z 2 ( r > fi ) are two tra- 
jectories which are close in the sense that |zj(r;s) — Z 2 (r;s)| < c(r)sl )- 7 for 
r G [0, 1]. Then 

|®(zf(r;s))-\F( 24 (r;s))| < [ \ £ s ii F i(4( r > *)) “ ^(4( r ! *)))l dr 

Jo 

< c(r|s|)-’ + 1 L/(> + 1 ), 

and so there is some fixed constant C such that 

W +1 . 


as desired. In particular, we have 

|z m (l,s)-$ F (s)|<C|s| i+1 . 

It remains to estimate |x J+1 (r;s) — z J+1 (r;s)|. Let 

|r°(r; s) - z°(r;s)| = 0, * J+1 (r;*) = *V( r\s )). 

Assuming that |z J (r;s) — r^(r; s) | < C(r|s|) J+1 , we estimate 
|i J+1 (r;s)-z J+1 (f;s)l = |$ i (r J '(r;s))-^( 2 J (r;s))| 

< |*V'(r;s))-*V(r; s »l 

+|*->'(z J (r; s)) - *(z J (r; s))| 

< J r I T j (%2 SiFi(z*(r\ s)) - T j 53 SiFi{z } (r- s)))| dr 

+ jT I T i (53«i F i( zi ( r ; *))“ 53s,fi(z J (r;s)))|dr 

< J' \ Tj C£ *)) ~ £ •* W(n *)))! dr 

+ ricE^ + vr 

Jo 

< r Ws\y + 2 

- C_ 7+2 _ ' 
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In particular, we have 


|x m (l,s) - 2 m (l,s)| < C|s| i+1 . 

Since $™(s) = x m (l;s), the lemma follows. ■ 

5 The main theorem 

In this section, we combine the lemmas of the two previous sections to prove 
that the easily computed approximating map A m maps trajectories of the 
explicitly integrable system £ to trajectories of the arbitrary system F with 
error 0(s m+1 ). 

Theorem 7 . Given a positive constant M and 6, a positive integer m i 
and a positive time T, there exists a positive constant C such that , given 
measurable controls tii(f) and ti 2 (t) satisfying |ti*(t)| < M for all t € [0 ,T], 
vector fields F\ and F 2 satisfying the BCE analyticity condition (1), and 
solutions x(t) and y(t) of 

i{t) = ui(t)Fi(T(t)) + u 2 (t)F 2 (x(t)), ar(0) = 0 

y(0 = u i + u 2 (t)E 2 (y(t)), y(0) = 0,< e [0.T], 

then |A m (y(<)) — z(t)| < ( Ct) m+l fort G [0,T], 

Proof. It is sufficient to show that the estimate holds at time T. Applying 
Lemma 5 twice, once to the system € and once to the system T shows that 

|x(T)-A(y(T))|<(aTr +1 . 

Also, we know from Lemma 6 that 

|x"‘(l;s)-x(l;s)|<C|sr +1 . 

Since A m (y(T)) = x m (l;s) and |s| < constant • T, the theorem follows. I 

6 Simultaneous Integration of Trajectories 

In this section, we describe an algorithm to integrate simultaneously a neigh- 
borhood of trajectories around a given fixed trajectory. Fix controls ui and 
U 2 and an initial condition x° £ F, and let z(t) denote the corresponding 
trajectory of the system £ . We give an algorithm yielding p trajectories in 
a tubular neighborhood around the fixed reference trajectory x(t). 
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Step 1. Solve the equation A(y) = x° for y. Let y° € T denote the root. 

Step 2. For points y 1 , . . y 9 in a neighborhood of y°, compute the trajec- 
tory y*{t) satisfying 

y(t) = + u 2 (t)E 2) y(0) = y*. 

Step 3. For each point y*(f) along the trajectory, compute the correspond- 
ing point A m (y>(*)) f for j = 1, . . p . 

The following theorem follows immediately from Theorem 7 . 

Theorem 8 For points y* sufficiently close to y° and small enough time 
t f the trajectories A m {y*{t)) all lie within a fixed tubular neighborhood of 
the reference trajectory x(t). These trajectories can be computed by simply 
evaluating the map A m , which can be precomputed , at different points . 


7 Examples 


We begin by showing why we require a condition such as BCH analyticity 
for the vector fields F,. 

Let ddd 


9X2 




where <^(xi,X2) is identically 1 if \x\ — Z2I > 1/4, and it is smooth and 
between 0 and 1 elsewhere, positive away from the diagonal, and vanishing 
to infinite order on the diagonal. Now the Z3 coordinate of exp(ay)oexp(aAf) 
is negative for all a > 0. But all brackets of X and Y vanish on the diagonal. 
Therefore, no formula of the form exp(u AT )oexp (ay) = exp(X+Y + brackets 
in X and Y can be valid. 

Figures 1 and 2 contain the Matbematica code we used to test the al- 
gorithm. Figure 3 contains the vector fields Ei and Fi and the A m map. 
Note that the map A m is a polynomial map. Figure 4 contains the result 
of flowing along the explicitly integrable flow in the C space, applying the 
map A m , and comparing the result to flowing in the T space using a Runge 
Kutta flow. 


8 



Brae [?. » wj : “Block [{i, j} t Table [Su*[w[[i]] D[w[[j]] ,x[i]] 
-w[[i]] D[T[[j]],x[i]] ,{i, Length W>] ,{j .Length [*]}] ] 

Flow [p_ • ▼_ * {rt _ , r0_ ,r 1 _}] s* Block[{xt t Tl ,xl ,i,j ,n>, 
n-Length[p] ; 
wl“v; 
xt*p; 

For [i“l , i<n , i++. 

xt[[i]] * Integrate [wl[[i]] ,rt]+p[[i]] ; 

Tl[[i+1]] * w[[i+l]]/.Table[x[j]->xt[[j]] ,{j»l,i>] 

Is 

xt [ [n] ] * Integrate [▼ 1 [[nj] ,rt]+p[DO] J 
(rt/ . {rt->rl})-(xt/ . {rt->rO» 

] 

Flow [p_ , ▼_ , {rt. »r 1_}] :“Flow[p,w,{rt ,0,ri}] 

Flow [p_ , w_ , {rt _}] :-Flow [p f ▼ , {rt , 0 , 1}] 

Flow[p_,wJ :“Flow[p,T ,{rt f 0, 1}] 


Figure 1: Mathematics code to compute brackets and integrable flows. 
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Picard[f_,m_] :«Block[{a,i, j ,r,n} t 

n»Length[f]; (* n * dimension *) 

a*f/. Table [x[i]->0,{i,nj] ; (* substitue x»0 *) 

For[j-l, j<*m, j++, (* Loop through ■ times... *) 

a* Integrate [a, r] ; 

a*f /. Table [x [i] - >a [ [i] ] ,{i,n>] ; (* substitue*) 

(* approximate by something integrable *) 

(* at increasing accuracy *) 

a-Table [formal [Series [a [[i]] ,{r,0, j-1}]] ,{i,n>] 

]; 

Integrate [a # {r,0,l}] (* Get the final answer *) 


Lambda [e_ , f _ , n_ , ■_] : -Block [{i , phiF , t , phiE , phiEInv} , 

C* e, f are functions whose values are vector fields *) 

(* n is the number of vf *s, m is the degree of approx. *) 
(* Get the F flow by Picard iteration *) 
phiF-Picard[Sum[t[i]f [iD ,{i,n}3.*3 ; 

(* Get the E flow by symbolic integration *) 
phiE-Flow [Table [0,{i,n}] ,Sum[t[i]e[i] ,{i,n}3] ; 

(* Invert the E flow symbolically *) 
phiEInv-Solve [Table [b [i] — phiE [ [i) 3 ,{i,n}] , 
Table[t[i] t {i,n>]][Cl]]; 

(♦ Form the composition *) 
phiF/.phiEInv 


Figure 2: Mathematica code to compute the map A m . 
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e[l]-{1.0,0} 
e[2]-{0.1,-x[l]} 
e [3] -Brae [e [2] ,e[l]] 
f [l]-{Sin[x[3]] ,0,x[2] Cos[x[3]]} 
f [2]-{0,Co*[x[3]] ,x[l]Sin[x[3]]} 
f [3] -Brae [f [2] ,f [1]] 

L3-Lanbda[e, 1,3,3] ; 

2 

b[l] b[2] b[l] (b[l] b[2] + 2 b[3]) 

{ + 

6 4 


2 

b[2] (b[l] b[2] + 2 b[3] ) 
b[2] + 


24 


b[l] b[2] b[l] b[2] + 2 b[3] 

+ — — 

2 2 


3 

(b[l] b[2] ♦ 2 b[3]) 
12 


Figure 3: The vector fields and the A m map. 

Flo v [{0 , 0 , 0} , (l+t+t * 2+t * 3) e [1] + (Sin [20t] ) e [2] , {t , . 1>] 

{0.105358, 0.0708073, -0.0045167} 

L3/ . Table [b[i] ->X [[i] ] , {i ,3}] 

{0.0000895593, 0.0708073, 0.00294345} 

RungeKuttaHiD[Join[{l}, (l+t+t*2+t"3)f [l]+(Sin{20t]) f[2]], 

{t,x[l] ,x{2] , i [3]}, {0,0, 0,0}, 0.02, 5] 

{0.1, 0.0000842708, 0.0708076, 0.00294349} 

Figure 4: Comparing the map \ m and a Runge Kutta flow. 
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